#!/usr/bin/env python
#-*- coding: utf-8 -*-

"""
Support for Avl (Athena Vortex Lattice)
See: http://web.mit.edu/drela/Public/web/avl/
"""
import os
import sys
import re
import pexpect
import logging
import numpy as np

import pat.utils as pu

patlogger = logging.getLogger('PAT')

def run_st(infile, outfile, case, debug=False):
    patlogger.info("Running st on %s runcase %d to %s", infile, case, outfile)
    if os.path.isfile(outfile): os.unlink(outfile)
    child = pexpect.spawn ('avl '+infile)
    if debug: child.logfile = sys.stdout
    child.expect ('AVL   c>')
    child.sendline ('OPER')
    child.expect ('OPER \(case \S+\)   c>')
    child.sendline ('{:d}'.format(case))
    child.expect ('OPER \(case \S+\)   c>')
    child.sendline ('X')
    child.expect ('OPER \(case \S+\)   c>')
    child.sendline ('ST')
    child.expect ('Enter filename, or <return> for screen output   s>')
    child.sendline (outfile)
    child.expect ('OPER \(case \S+\)   c>')
    if debug: print str(child)


def run_eig(infile, outfile, case, debug=False):
    if os.path.isfile(outfile): os.unlink(outfile)
    child = pexpect.spawn ('avl '+infile)
    if debug: child.logfile = sys.stdout
    child.expect ('AVL   c>')
    child.sendline ('MODE')
    child.expect ('MODE   c>')
    child.sendline ('{:d}'.format(case))
    child.expect ('MODE   c>')
    child.sendline ('S')
    child.expect ('Enter output filename \(or <Return>\):')
    child.sendline (outfile)
    child.expect ('MODE   c>')
    if debug: print str(child)


def parse_st(filename, debug=False, surfaces=('elevator', 'aileron')):
    """
    Extracts coefficients from an avl st file
    """
    parser = [(r'.*Sref =\s+(\S+)\s+Cref =\s+(\S+)\s+Bref =\s+(\S+)',('Sref', 'Cref', 'Bref')),
              (r'.*Alpha =\s+(\S+)\s+', ('Alpha',)),
              (r'.*CYtot =\s+(\S+)\s+ Cmtot =\s+(\S+)\s+',('CYtot', 'Cmtot')),
              (r'.*CLtot =\s+(\S+)\s+',('CLtot',)),
              (r'.*CDtot =\s+(\S+)\s+',('CDtot',)),
              (r'.*CDvis =\s+(\S+)\s+ CDind =\s+(\S+)\s+',('CDvis', 'CDind')),
              (r'.*CLff  = \s+(\S+)\s+ CDff  =\s+(\S+)\s+',('CLff', 'CDff')),
              # surfaces will be inserted here...              
              (r'.*CLa =\s+(\S+)\s+ CLb =\s+(\S+)\s+', ('CLa', 'CLb')),
              (r'.*CYa =\s+(\S+)\s+ CYb =\s+(\S+)\s+', ('CYa', 'CYb')),
              (r'.*Cla =\s+(\S+)\s+ Clb =\s+(\S+)\s+', ('Cla', 'Clb')),
              (r'.*Cma =\s+(\S+)\s+ Cmb =\s+(\S+)\s+', ('Cma', 'Cmb')),
              (r'.*Cna =\s+(\S+)\s+ Cnb =\s+(\S+)\s+', ('Cna', 'Cnb')),
              (r'.*CLp =\s+(\S+)\s+ CLq =\s+(\S+)\s+ CLr =\s+(\S+)\s+', ('CLp', 'CLq', 'CLr')),
              (r'.*CYp =\s+(\S+)\s+ CYq =\s+(\S+)\s+ CYr =\s+(\S+)\s+', ('CYp', 'CYq', 'CYr')),
              (r'.*Clp =\s+(\S+)\s+ Clq =\s+(\S+)\s+ Clr =\s+(\S+)\s+', ('Clp', 'Clq', 'Clr')),
              (r'.*Cmp =\s+(\S+)\s+ Cmq =\s+(\S+)\s+ Cmr =\s+(\S+)\s+', ('Cmp', 'Cmq', 'Cmr')),
              (r'.*Cnp =\s+(\S+)\s+ Cnq =\s+(\S+)\s+ Cnr =\s+(\S+)\s+', ('Cnp', 'Cnq', 'Cnr')),
              # coef for surfaces will be appended here
              ]
    for sfc in reversed(surfaces):
        parser.insert(7, (r'.*{:s}\s+ =\s+(\S+)\s+'.format(sfc), (sfc,)))

    for coef in ['CL', 'CY', 'Cl', 'Cm', 'Cn', 'CDff', 'e']:
        regexp = ".+"+"".join(['{:s}d{:d} =\s+(\S+)\s+'.format(coef, i+1) for i,s in enumerate(surfaces)])
        var = tuple(['{:s}d{:d}'.format(coef, i+1) for i,s in enumerate(surfaces)])
        parser.append((regexp, var))
    
    f = open(filename)
    i = 0; values = {}
    for line in f:
        patlogger.debug("vehicles.fixed_wing.avl.parse_st: %s", line.replace('\n', ''))
        if i < len(parser):
            m=re.match(parser[i][0], line)
            if m:
                patlogger.debug("vehicles.fixed_wing.avl.parse_st: match found %s %s", parser[i][1], m.groups())
                for j, v in enumerate(m.groups()):
                    values[parser[i][1][j]] = v
                i += 1
    patlogger.info("vehicles.fixed_wing.avl.parse_st: matched %d/%d regular expressions in %s",
                   i, len(parser), filename)
    return values


def get_masses(filename, debug=False):
    """
    Runs Avl for parsing a mass file.
    returns computed mass and inertia
    """
    cfg_dir = os.path.dirname(filename) # avl needs to run from the config files directory
    cur_dir = os.getcwd()               # save that to be able to return
    os.chdir(cfg_dir)
    cfg_file = os.path.basename(filename)
    child = pexpect.spawn ('avl '+cfg_file)
    if debug: child.logfile = sys.stdout
    child.expect(r'Mass \s+=\s+(\S+)\s+kg')
    (mass,) = child.match.groups()
    child.expect(r'Ixx -Ixy -Ixz   \|\s+(\S+)\s+(\S+)\s+(\S+)')
    (Ixx, mIxy, mIxz) = child.match.groups()
    child.expect(r'Iyy -Iyz = \|\s+(\S+)\s+(\S+)')
    (Iyy, mIyz) = child.match.groups()
    child.expect(r'Izz   \|\s+(\S+)')
    (Izz,) = child.match.groups()
    child.expect ('AVL   c>')
    child.sendline ('Quit')
    child.expect(pexpect.EOF)
    I = np.array([[ float(Ixx),  float(mIxy), float(mIxz)],
                  [ float(mIxy), float(Iyy),  float(mIyz)],
                  [ float(mIxz), float(mIyz), float(Izz)]])
    res = {'mass':float(mass), 'inertia':I}
    os.chdir(cur_dir)
    return res



################################################################################
#
#  Below, crap...
#
#








def parse_mass(filename):
    """
    Extracts coefficients from an avl mass file - warning, totally uncomplete
    """
    res = {'filename':filename, 'masses':[]}
    f = open(filename)
    #re_num = r'[-+]?\d*\.\d+|\d'
    re_num = r'[-+]?\d*\.\d+[E[-]?\d+]?|\d'
    regex_cmt = r'^.*#.*'
    regex_mass = r'\s*('+re_num+')\s+('+re_num+')\s+('+re_num+')\s+('+re_num+')\s+('+re_num+')\s+('+re_num+')\s+('+re_num+')\s+.*'
    for line in f:
        if not re.match(regex_cmt, line):
            match=re.match(regex_mass, line)
            if match:
                patlogger.info("vehicles.fixed_wing.avl.parse_st: match found %s", match.groups())
                (m, x, y, z, ix, iy, iz) = match.groups()
                res['masses'].append((float(m), (x, y, z), (float(ix), float(iy), float(iz))))
    res['mass'] = sum([m for (m, l, i) in res['masses']])
    res['inertia'] = res['masses'][0][2]
    return res



def print_xml(outfile, ac_name, values_CL0, values_CD0, values_der):
    f = open(outfile, 'w')
    f.write('''<?xml version="1.0"?>
<fdm_config name="{:s}">
  <masses>
    <mass unit="kg"> 0.6 </mass>
    <inertia unit="kg.m2">
0.010, 0.,    0.
0.,    0.022, 0.
0.,    0.,    0.025
    </inertia>
  </masses>
 
  <propulsion>
    <engine name="engine" pos="0., 0., 0." align="0., 0." fmax="4" rhoi="1.225" nrho="1." Vi="15." nV="0." tau="0.1"/>
  </propulsion>
'''.format(ac_name))
    f.write(r'''
  <aerodynamics>
    <Sref unit="m2" comment="Reference surface">{:f}</Sref>
    <Cref unit="m"  comment="Reference chord">  {:f}</Cref>
    <Bref unit="m"  comment="Reference length"> {:f}</Bref>
'''.format(float(values_der['Sref']), float(values_der['Cref']), float(values_der['Bref'])))
    f.write(r'''    <stability>
       <misc CL0="{:f}" alpha0="{:f}" CD0="{:f}" CD_k1="{:f}" CD_k2="{:f}" Cm0="{:f}"/>
'''.format(float(values_CL0['CLtot']),
           0., #pu.rad_of_deg(float(values_der['Alpha'])),
           float(values_CD0['CDtot']),
           0.,
           (float(values_der['CDtot'])-float(values_CD0['CDtot']))/float(values_der['CLtot'])**2,
            float(values_CL0['Cmtot'])))

    f.write(\
r'''       <alpha                   CL="{: f}" CY="{: f}"                Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(float(values_der['CLa']), float(values_der['CYa']),
           float(values_der['Cla']), float(values_der['Cma']), float(values_der['Cna'])))
    f.write(\
r'''       <beta                    CL="{: f}" CY="{: f}"                Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(float(values_der['CLb']), float(values_der['CYb']),
           float(values_der['Clb']), float(values_der['Cmb']), float(values_der['Cnb'])))
    f.write(\
r'''       <p                       CL="{: f}" CY="{: f}"                Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(float(values_der['CLp']), float(values_der['CYp']),
           float(values_der['Clp']), float(values_der['Cmp']), float(values_der['Cnp'])))
    f.write(\
r'''       <q                       CL="{: f}" CY="{: f}"                Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(float(values_der['CLq']), float(values_der['CYq']),
           float(values_der['Clq']), float(values_der['Cmq']), float(values_der['Cnq'])))
    f.write(\
r'''       <r                       CL="{: f}" CY="{: f}"                Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(float(values_der['CLr']), float(values_der['CYr']),
           float(values_der['Clr']), float(values_der['Cmr']), float(values_der['Cnr'])))
    f.write(\
r'''       <surface name="aileron"  CL="{: f}" CY="{: f}" CD="{: f}" Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(pu.deg_of_rad(float(values_der['CLd2'])), pu.deg_of_rad(float(values_der['CYd2'])), pu.deg_of_rad(float(values_der['CDffd2'])), 
           pu.deg_of_rad(float(values_der['Cld2'])), pu.deg_of_rad(float(values_der['Cmd2'])), pu.deg_of_rad(float(values_der['Cnd2']))))
    f.write(\
r'''       <surface name="elevator" CL="{: f}" CY="{: f}" CD="{: f}" Cl="{: f}" Cm="{: f}" Cn="{: f}"/>
'''.format(pu.deg_of_rad(float(values_der['CLd1'])), pu.deg_of_rad(float(values_der['CYd1'])), pu.deg_of_rad(float(values_der['CDffd1'])), 
           pu.deg_of_rad(float(values_der['Cld1'])), pu.deg_of_rad(float(values_der['Cmd1'])), pu.deg_of_rad(float(values_der['Cnd1']))))
            
    f.write(r'''    </stability>
  </aerodynamics>
</fdm_config>
''')
    f.close()




#
#
#
#m is to modify the parameters first
#Murat • Now
#ok
#and then C1 and C2 ?
#Now
#there are some weird stuff when you change some parameters, like if you change the mass in M
#and forget to go to C1 for velocity

#but then when you run C1 or C2 it overwrites the parameters
#yes like gravity, mass, etc should have to be done in M
#thats why if you change the mass, you have to re enter C1 and re write the V
#so that it calculates the new CL coeff
#
#  K2 = 1/(pi*Bref**2/Sref*e)
#

#and i still need to figure which coefficients are degrees 😃
#:) hehe :)
#yes that is a little bit stupid on AVL output
#on an stability derivative output, which you obtain with the ST command in AVL oper mode, only the control surfaces are given in degrees, all of the other stuff is in radians
#Murat • 14 mins
#so Clq
#q is rad/s
#and Cla
#alpha is radians?
#yes
#Murat • 13 mins
#but alpha0 is degres....
#CLd1 is in degrees
#Murat • 13 mins
#Alpha =   4.50979
#i assume this degres
#and yest Cla is radians?
#12 mins
#
#
#yes that is before the derivative part on the output
#Murat • 12 mins
#ok.... I'm adjusting that
#11 mins
#
#
#yes Cla is in radians
#Murat • 12 mins
#it's a factor of 60... we should be able to see if the response is off by that much 😃
#11 mins
#
#
#and elevator and aileron values for the equilibrium is in degrees, looks like all the informative stuff is given in degrees just before the derivative part.
#You should be also careful for the deflections, on the avl file you can define a coefficient which is multiplied by the control input.
#Murat • 7 mins
#mmmmm
#6 mins
#
#
#I am usually taking it as 0.75 in the avl file. Which means that when you give 10 degrees of control deflection input to the model, the aerodynamic model only moves 7.5degrees and gives the forces and moments according to that 7.5degrees
#Murat • 5 mins
#why would you do that?
#3 mins
#
#
#This is because of the simplified model used in AVL for the control surfaces. It uses the camber change, and generates a new camber line for the wing when you give a control surface deflection. As the direct camber change in the AVL model is giving a little bit optimistic force and moment generation compared to real life, we have to reduce it to make it close to reality...
#In real life there will be more flow separation which results less control surface efficiency, so it has to be taken into account...
#But it really depends on the manufacturing of the wing, the airfoil's pressure loading along the control surface, etc… and I think there is no "real" way of estimating it correctly, at least not that I am aware of.
#
#    I ahem seen that for foam airplanes 0.70-0.75 looks ok, and for composite sharp trailing edged models maybe something like 0.80-0.85 would be better
#0.75 number was not really experimental or experience in the beginning, I think it was in the schligting book for boundary layer theory, where he was suggesting this efficiency factor, so if we find the statement in the book, you can use it as a reference :)
    

#for drag we write: CD = CD0 + k1.CL + K2.CL^2 + Ksfc.delta_sfc
#CD0 is what is called CDtot in the AVL file
#what are K1 and K2? CDvis, CDind, CDff ?
#15 mins


#CDtot of AVL is the total drag, including everything normally. But you have to be aware of that AVL doesn't take into account the viscous drag. so CDvis has to be added correctly to the AVL output. In order to do that you either have to know the exact working conditions that you are going to do the analysis (flight speed, air density, kinematic viscosity of air) so that you can define a polynomial curve for for your airfoils and AVL uses these for the CDvis calculations. Or you just define an approximative value for the CDp on the airplane .avl file which will remain constant for every condition...
#The latter can be acceptable if you are only doing some conceptual design trials
#Murat • 7 mins
#haha 😃
#6 mins


#but if you really need something precise (actually nothing is going to be that much precise :) )then you have to work on the airfoils and define the polynomial fit curve coefficients...
#Murat • 6 mins
#so..... for something not precise... CD0 = CDtot and K1=0 and k2=0 ?
#5 mins


#For something not precise, you can run the AVL for zero lift condition, in oper , you have to set the elevator controls to zero pitching moment by writing "D2 PM 0" assuming that the elevator control surface is D2. and then still in oper mode, you have to set the alpha value to the CL=0 by writing "a c 0", then run the condition by "x" and look at the results on the trefftz plane by "t" command.

#there CDtot which is written CD on the treftz plane plot, is the CD0,

#then for K1 and K2 values… :
#well there is nothing easy for the first order K1 value, which is usually not significant...
#but for K2 you can use the oswald efficiency factor which is "e"
#the formula is (CL^2)/(pi*AR*e)
#for the K2*CL^2

#but..... we have an alpha0
#i mean we write
#CL = CL0 + CLa*(alpha-alpah0)+...
#and CD = CD0 + K1.CL + K2.CL^2 +...
#so in this case isn't CL0 what is called CLtot and CD0 what is called CDtot ?
#and alpha0 what is called ALpha in the avl file?
#3 mins


#so exactly, when we set the alpha to zero CL angle with "a c 0" K1 and K2 will be cancelled out to find the CD0 which is equal to CDtot for that alpha. For the K2 you have to run another alpha which you are mostly interested in (such as cruise aplha) and get the e value for the K2 coefficient calculation from that new run.
#Murat • 1 min

#i think I understand? I will try to write that
#sorry I have to explain it in more detail,
#AVL is able to calculate almost everything(except the viscous effects) so it doesn't need any K1 K2 or CD0 values...
#
#But in analytic calculations, we usually have to simplify it to the Cd= Cd0+k1Cl+k2Cl^2+...
#AVL calculates everything in one, and gives the Cd = Cdtot
#the same for CL= CLtot
#Murat • 1 min
#ok... so if I make the two points identification you suggest, later i can run more avl cases and they should more or less match what my model gives ?
#Now


#in analytic CL= CL0+CLa*a+CLelevator*elevator+…..

#
# yes AVL is linear for most of the values, so 2 case runs will be sufficient for a simple model.
# normally 1 case is sufficient for the derivatives
#
# then 2nd case will be used for CD0
#yes! it starts to make sense
#thank you! back to work 😃
#Now


#hmm maybe a 3rd case for Alpha zero and CL0
#Murat • Now

#what? 😃
#CL0 is the CL value for angle of attck=0 so this needs 1 run
#CD0 is the CDtot when CL=0 so this is another alpha, so this is 2nd run
#
#and then 3rd run when alpha corresponds to CL=CLcruise, this is for the derivatives
#Murat • Now
#on these 3 cases ptcihng moment has to be constrained by the elevator I think
#ok!
